rm(list=ls())
graphics.off()

library(readxl)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(ggpubr)
library(rstatix)
library(nlme)
library(matrixcalc)
library(Biodem)
library(stats)
library(multcomp)
library(hrbrthemes)
library(viridis)
library(gapminder)
library(devtools)
library(stringr)
library(patchwork)
library(gridExtra)
library(cowplot)
library(extrafont)
library(gtable)
library(ggtext)
library(effsize)
library(data.table)
library(psych)
library(psychotools)
library(psychTools)
library(BlandAltmanLeh)

######################### START ###########################################

setwd("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3")

SR <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/SucessRateMean.csv")

Rep <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/PercentageMean.csv")


colnames(SR) <- c("ID","aclow","bcchall","cchigh","dmlow","emchall","fmhigh")
colnames(Rep) <- c("ID","aclow","bcchall","cchigh","dmlow","emchall","fmhigh")

###SR 
SR_long <- pivot_longer(SR, cols=c("aclow","bcchall","cchigh","dmlow","emchall","fmhigh"))
div <- rep(c("Mental","Motor"), each = 3)
SR_long$div <- rep(div,nrow(SR)) 

SR_long = split(SR_long, SR_long$div)

SR_coglong = SR_long$Mental[,1:3]

SR_motlong =  SR_long$Motor[,1:3]

###

plot <- ggplot(SR_long, aes(name, value, ID, fill = name)) +
  geom_boxplot() +
  labs(x = 'Intensity Level',
       title = 'SR',
       y = "SR") +
  theme(plot.title = element_text(hjust = 0.5, size = 20),
        legend.title=element_blank(),
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "grey", size = 0.5),
        axis.title = element_text(size = 13),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 13)) +
  scale_x_discrete(label = element_blank()) +
  scale_y_continuous(limits = c(-0.015,0.223)) +
  guides(fill="none") +
  scale_fill_manual(values=c("steelblue4", "lightblue1", "cadetblue3", "steelblue4", "lightblue1", "cadetblue3")) +
  annotate(geom = "text", x = 0.83, y = -0.009, label = "Very Low", hjust = 0, vjust = 0, size = 3, color = "gray50") +  
  annotate(geom = "text", x = 0.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.009, label = "Challenging", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.85, y = -0.009, label = "Very High", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.87, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  facet_grid(. ~ div, scales = "free", space = "free")


#ggsave("HRVplot.png", plot, width = 10.5, height = 6.92)

###Statistics

###Test Normality
Norm.SR.CogEasy = shapiro.test(SR$aclow)
Norm.SR.CogMiddle = shapiro.test(SR$bcchall)
Norm.SR.CogDifficult = shapiro.test(SR$cchigh)

Norm.SR.MotEasy = shapiro.test(SR$dmlow)
Norm.SR.MotMiddle = shapiro.test(SR$emchall)
Norm.SR.MotDifficult = shapiro.test(SR$fmhigh)

###Perform Test and Post-Hoc
### Cognitive 
if (Norm.SR.CogEasy$p.value > 0.05 & Norm.SR.CogMiddle$p.value > 0.05 & Norm.SR.CogDifficult$p.value > 0.05){
  print(anova_test(data=SR_coglong, dv = value, wid=ID, within = name))
  pairwise.t.test(SR_coglong$value, SR_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=SR_coglong, value ~ name|ID))
  pairwise.wilcox.test(SR_coglong$value, SR_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}

###Motor
### Cognitive 
if (Norm.SR.MotEasy$p.value > 0.05 & Norm.SR.MotMiddle$p.value > 0.05 & Norm.SR.MotDifficult$p.value > 0.05){
  print(anova_test(data=SR_motlong, dv = value, wid=ID, within = name))
  pairwise.t.test(SR_motlong$value, SR_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=SR_motlong, value ~ name|ID))
  pairwise.wilcox.test(SR_motlong$value, SR_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}


#ICC analysis

data_SR_Test <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/SucessRateDay1.csv")
data_SR_ReTest <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/SucessRateDay2.csv")


## SR

## Cognitive

### Easy
data_SR_Cog_easy = data.table(data_SR_Test$Cog_Easy, data_SR_ReTest$Cog_Easy)
#data_SR_easy_RRn = data_SR_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SR_Cog_easy =ICC(data_SR_Cog_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SR_Cog_easy<- bland.altman.plot(data_SR_Test$Cog_Easy,data_SR_ReTest$Cog_Easy, main="Test-Retest SR Easy", xlab="Means", ylab="differences")

print(ICC_SR_Cog_easy)

###Middle
data_SR_Cog_middle = data.table(data_SR_Test$Cog_Optimal, data_SR_ReTest$Cog_Optimal)
#data_SR_Cog_middle = data_SR_Cog_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SR_Cog_middle =ICC(data_SR_Cog_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SR_Cog_middle <-bland.altman.plot(data_SR_Test$Cog_Optimal,data_SR_ReTest$Cog_Optimal, main="Test-Retest SR Optimal", xlab="Means", ylab="differences")

print(ICC_SR_Cog_middle)


###Difficult
data_SR_Cog_difficult = data.table(data_SR_Test$Cog_Difficult, data_SR_ReTest$Cog_Difficult)
#data_SR_Cog_difficult = data_SR_Cog_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SR_Cog_difficult =ICC(data_SR_Cog_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SR_Cog_difficult<- bland.altman.plot(data_SR_Test$Cog_Difficult,data_SR_ReTest$Cog_Difficult, main="Test-Retest SR Difficult", xlab="Means", ylab="differences")

print(ICC_SR_Cog_difficult)

## Motor

### Easy
data_SR_Mot_easy = data.table(data_SR_Test$Mot_Easy, data_SR_ReTest$Mot_Easy)
#data_SR_easy_RRn = data_SR_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SR_Mot_easy =ICC(data_SR_Mot_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SR_Mot_easy<- bland.altman.plot(data_SR_Test$Mot_Easy,data_SR_ReTest$Mot_Easy, main="Test-Retest SR Easy", xlab="Means", ylab="differences")

print(ICC_SR_Mot_easy)

###Middle
data_SR_Mot_middle = data.table(data_SR_Test$Mot_Optimal, data_SR_ReTest$Mot_Optimal)
#data_SR_Mot_middle = data_SR_Mot_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SR_Mot_middle =ICC(data_SR_Mot_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SR_Mot_middle <-bland.altman.plot(data_SR_Test$Mot_Optimal,data_SR_ReTest$Mot_Optimal, main="Test-Retest SR Optimal", xlab="Means", ylab="differences")

print(ICC_SR_Mot_middle)


###Difficult
data_SR_Mot_difficult = data.table(data_SR_Test$Mot_Difficult, data_SR_ReTest$Mot_Difficult)
#data_SR_Mot_difficult = data_SR_Mot_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_SR_Mot_difficult =ICC(data_SR_Mot_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_SR_Mot_difficult<- bland.altman.plot(data_SR_Test$Mot_Difficult,data_SR_ReTest$Mot_Difficult, main="Test-Retest SR Difficult", xlab="Means", ylab="differences")

print(ICC_SR_Mot_difficult)







###Rep 
Rep_long <- pivot_longer(Rep, cols=c("aclow","bcchall","cchigh","dmlow","emchall","fmhigh"))
div <- rep(c("Mental","Motor"), each = 3)
Rep_long$div <- rep(div,nrow(Rep)) 

Rep_long = split(Rep_long, Rep_long$div)

Rep_coglong = Rep_long$Mental[,1:3]

Rep_motlong =  Rep_long$Motor[,1:3]

###

plot <- ggplot(Rep_long, aes(name, value, ID, fill = name)) +
  geom_boxplot() +
  labs(x = 'Intensity Level',
       title = 'Rep',
       y = "Rep") +
  theme(plot.title = element_text(hjust = 0.5, size = 20),
        legend.title=element_blank(),
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "grey", size = 0.5),
        axis.title = element_text(size = 13),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 13)) +
  scale_x_discrete(label = element_blank()) +
  scale_y_continuous(limits = c(-0.015,0.223)) +
  guides(fill="none") +
  scale_fill_manual(values=c("steelblue4", "lightblue1", "cadetblue3", "steelblue4", "lightblue1", "cadetblue3")) +
  annotate(geom = "text", x = 0.83, y = -0.009, label = "Very Low", hjust = 0, vjust = 0, size = 3, color = "gray50") +  
  annotate(geom = "text", x = 0.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.009, label = "Challenging", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.85, y = -0.009, label = "Very High", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.87, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  facet_grid(. ~ div, scales = "free", space = "free")


#ggsave("HRVplot.png", plot, width = 10.5, height = 6.92)

###Statistics

###Test Normality
Norm.Rep.CogEasy = shapiro.test(Rep$aclow)
Norm.Rep.CogMiddle = shapiro.test(Rep$bcchall)
Norm.Rep.CogDifficult = shapiro.test(Rep$cchigh)

Norm.Rep.MotEasy = shapiro.test(Rep$dmlow)
Norm.Rep.MotMiddle = shapiro.test(Rep$emchall)
Norm.Rep.MotDifficult = shapiro.test(Rep$fmhigh)

###Perform Test and Post-Hoc
### Cognitive 
if (Norm.Rep.CogEasy$p.value > 0.05 & Norm.Rep.CogMiddle$p.value > 0.05 & Norm.Rep.CogDifficult$p.value > 0.05){
  print(anova_test(data=Rep_coglong, dv = value, wid=ID, within = name))
  pairwise.t.test(Rep_coglong$value, Rep_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=Rep_coglong, value ~ name|ID))
  pairwise.wilcox.test(Rep_coglong$value, Rep_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}

###Motor
### Cognitive 
if (Norm.Rep.MotEasy$p.value > 0.05 & Norm.Rep.MotMiddle$p.value > 0.05 & Norm.Rep.MotDifficult$p.value > 0.05){
  print(anova_test(data=Rep_motlong, dv = value, wid=ID, within = name))
  pairwise.t.test(Rep_motlong$value, Rep_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=Rep_motlong, value ~ name|ID))
  pairwise.wilcox.test(Rep_motlong$value, Rep_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}


#ICC analysis

data_Rep_Test <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/PercentageDay1.csv")
data_Rep_ReTest <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/PercentageDay2.csv")


## Rep

## Cognitive

### Easy
data_Rep_Cog_easy = data.table(data_Rep_Test$Cog_Easy, data_Rep_ReTest$Cog_Easy)
#data_Rep_easy_RRn = data_Rep_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_Rep_Cog_easy =ICC(data_Rep_Cog_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_Rep_Cog_easy<- bland.altman.plot(data_Rep_Test$Cog_Easy,data_Rep_ReTest$Cog_Easy, main="Test-Retest Rep Easy", xlab="Means", ylab="differences")

print(ICC_Rep_Cog_easy)

###Middle
data_Rep_Cog_middle = data.table(data_Rep_Test$Cog_Optimal, data_Rep_ReTest$Cog_Optimal)
#data_Rep_Cog_middle = data_Rep_Cog_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_Rep_Cog_middle =ICC(data_Rep_Cog_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_Rep_Cog_middle <-bland.altman.plot(data_Rep_Test$Cog_Optimal,data_Rep_ReTest$Cog_Optimal, main="Test-Retest Rep Optimal", xlab="Means", ylab="differences")

print(ICC_Rep_Cog_middle)


###Difficult
data_Rep_Cog_difficult = data.table(data_Rep_Test$Cog_Difficult, data_Rep_ReTest$Cog_Difficult)
#data_Rep_Cog_difficult = data_Rep_Cog_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_Rep_Cog_difficult =ICC(data_Rep_Cog_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_Rep_Cog_difficult<- bland.altman.plot(data_Rep_Test$Cog_Difficult,data_Rep_ReTest$Cog_Difficult, main="Test-Retest Rep Difficult", xlab="Means", ylab="differences")

print(ICC_Rep_Cog_difficult)

## Motor

### Easy
data_Rep_Mot_easy = data.table(data_Rep_Test$Mot_Easy, data_Rep_ReTest$Mot_Easy)
#data_Rep_easy_RRn = data_Rep_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_Rep_Mot_easy =ICC(data_Rep_Mot_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_Rep_Mot_easy<- bland.altman.plot(data_Rep_Test$Mot_Easy,data_Rep_ReTest$Mot_Easy, main="Test-Retest Rep Easy", xlab="Means", ylab="differences")

print(ICC_Rep_Mot_easy)

###Middle
data_Rep_Mot_middle = data.table(data_Rep_Test$Mot_Optimal, data_Rep_ReTest$Mot_Optimal)
#data_Rep_Mot_middle = data_Rep_Mot_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_Rep_Mot_middle =ICC(data_Rep_Mot_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_Rep_Mot_middle <-bland.altman.plot(data_Rep_Test$Mot_Optimal,data_Rep_ReTest$Mot_Optimal, main="Test-Retest Rep Optimal", xlab="Means", ylab="differences")

print(ICC_Rep_Mot_middle)


###Difficult
data_Rep_Mot_difficult = data.table(data_Rep_Test$Mot_Difficult, data_Rep_ReTest$Mot_Difficult)
#data_Rep_Mot_difficult = data_Rep_Mot_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_Rep_Mot_difficult =ICC(data_Rep_Mot_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_Rep_Mot_difficult<- bland.altman.plot(data_Rep_Test$Mot_Difficult,data_Rep_ReTest$Mot_Difficult, main="Test-Retest Rep Difficult", xlab="Means", ylab="differences")

print(ICC_Rep_Mot_difficult)

